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Abstract 

Expressions for intermolecular forces and torques, derived from pair potentials between rigid non- 
spherical units, are presented. The aim is to give compact and clear expressions, which are easily 
generalised, and which minimise the risk of error in writing molecular dynamics simulation pro- 
grams. It is anticipated that these expressions will be useful in the simulation of liquid crystalline 
systems, and in coarse-grained modelling of macromolecules. 
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I. INTRODUCTION 



This paper is intended to clarify the way in which forces and torques may be computed 
from an expression for the pair potential between two, generally nonlinear, molecules. These 
forces and torques are then typically incorporated into molecular dynamics simulation pro- 
grams The particular focus is on pair potentials which are expressed directly 
as functions of the molecular centre-centre vector and orientational variables such as Eu- 
ler angles, quaternion (Euler) parameters, or rotation matrix elements. Some examples, 
based on potentials used in the simulation of liquid crystals and macromolecules, are given, 
but the intention is also to provide a formalism that is easily generalised. The more com- 
mon case of atomic site-site potentials is not considered here. Several schemes exist for 
systematically "coarse-graining" such atomistically-detailed interaction potentials, often by 
performing simulations over relatively short times and fitting either coarse-grained pair dis- 
tribution functions or aggregated forces 

mm 

, but this is also not the focus of the 
present paper. The aim is to streamline, and make more robust, the process of analytically 
deriving forces and torques, once the coarse-grained potentials have been determined, and 
parametrised in a standard way. Even with the help of symbolic algebra packages, per- 
forming this task is prone to error and may result in molecular dynamics code which is 
inefficient, hard to understand, and difficult to adapt. By the use of some examples, based 
on spheroidal geometry, it is shown here that standard vector and matrix notation helps to 
simplify the formulae, leading to an improved insight into the machinery of the potentials 
as well as giving a more reliable route to the simulation code. 

Consider two nonlinear molecules A and B, position vectors r^, r^. Define the inter- 
molecular vector r = tab = fA — ^b, distance r = \r\, and unit vector r = r/r. Specify the 
orientations by the orthogonal rotation matrices a, b, which convert from space-fixed (xyz) 
to molecule-fixed (123) coordinates |9|: 



0-2 



^ 0,lx O^ly 0,1. ^ 
0'2x 0,2y 0.2, 

\o,3x a^y 0,32 y 



bix biy bi 

b2x b2y b2z 



\b3x bsy b^z J 



The rows of these matrices are the three molecule-fixed orthonormal principal axis vectors 



m = 1,2, 3. Note that it is possible to write [a]m/x 
without ambiguity. 



a 



mj n 



"mill 



fx 



X, y, z. 



It is assumed that the pair potential may be written as a function U = U{r,a,b) and 
that it is invariant to a global rotation of the coordinates. Quite generally the force on 
A due to B will be the negative of the force on B due to A, so 

f --f -f- ^ 

J A- JB = J--q^- 

Depending on the form of the potential it may be possible to evaluate this derivative directly, 
or it may be more convenient to separate the components along, and perpendicular to, the 
line of centres: 

^_ dUdr dU fdr\__dU^ ^_^dU ^ ^ 

Or Or Or \dr J dr dr 

On the right, constructs a dyadic matrix from two vectors, and 1 is the unit 3x3 matrix. 
The torque calculation follows Price et al. Q], with minor variations. Consider the negative 
of the derivative of U with respect to rotation of molecule A through an angle about any 
axis represented by a unit vector tp. By definition, this gives the corresponding component 
of the torque ta, and the chain rule may be used to express this in terms of the body-fixed 
unit vectors 

^ dU -syy dU dcLm dU ^ ^ ^ 

dip ^ dCira dip ^ dCLra ^ 

The last equation follows from the rotation formula 9] for any vector de/dip = x e. Cyclic 
permutation of the vectors in the scalar triple product then gives 

dU 



tl) ■ Ta = -'ll^ -^CLr, 



dCLr. 



Choosing ip to be each of the coordinate directions in turn allows the identification of every 
component of the torque, and the expression for is similar 

rA = -l^amX —- , TB = -}_^K,x . (2) 

In the succeeding sections, more specific examples of the above formulae will be given for 
particular choices of pair potential U . 
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II. SCALAR PRODUCT POTENTIALS 



A very common case, and a useful preamble for what follows, is when the pair potential 
between A and B may be written in the form ^| 

U = U{r, a, b) = U{r, {a„ ■ f}, {K ■ r}, {a^ • b„}) . (3) 

This is a function of the centre-centre separation r, and all possible scalar products of the 
unit vectors r, and b„. It is manifestly invariant to a global rotation of coordinates. The 
force may be written 

^ _ _dU_ _ _dl^dr_ _sp( dU d{a^ ■ r) ^ dU djbm ■ r) 
dr dr dr ^ \ d{drn ■ r) dr d{h„, ■ r) dr 

The sum ranges over all the orientation vectors on both molecules, {a,m, bm}- The notation 

= V — (v ■ r)r = —r x (f x 

is short for the component of a vector v perpendicular to r. The derivatives of the potential 
appearing in eqn are easily evaluated, assuming that it has the general form of eqn Q. 
The torque calculation once more follows Price et al. |ilO,]. The derivative of U with respect 
to rotation of molecule A through an angle ip about an axis '0 takes the form 

^. dU ^ dU d{r-am) ^ dU d{b 

All combinations of unit vectors for which one, a^, rotates with the molecule while the 
other, r or remains stationary, appear on the right. The rotation formula in this case 
gives d{v ■ am)/dil) = v - ip x = —'0 ■ v x and once again -0 may be chosen arbitrarily. 
The explicit result for molecules A and B is 



Note that the combination of eqns (jlj), (j5a|) and (j5b|) gives 

TA + TB + rx f = 0. (6) 



This is the expression of local angular momentum conservation, which follows directly from 
the invariance of the potential energy, eqn Q, to a rotation of all the vectors f, and 6„ 
together. 



III. BIAXIAL GAY-BERNE, AND RELATED, POTENTIALS 



The original Gay-Berne family of potentials 
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3,Q 



was developed for molecules 



having the symmetry, and approximate shape, of uniaxial ellipsoids. They have now been 
generalised to include biaxiality Q, Q 0, Q . The general form is \li 

Ugb = £i (a, b) 4(a, b, f) f/shapc {r, (r{r, a, b)) 

Here the two e terms represent orientation-dependent strength factors, raised to powers 
fi, V which characterise different variants of the model. The last term roughly defines the 
molecular shape, via an orient at ion- dependent distance function a which will be discussed 
in more detail below, f/shapc typically is based on a Lennard- Jones 12-6 potential 



shape ' ) 



r, cr 



4en 



12 



(7\6 

r 



(7) 



although, as shall be seen, this has been refined somewhat. The Gay -Berne potentials have 
been extensi vely used to model small organic molecules 3, 2^ 21 1 and liquid crystalline 



systems |2^ |23|; linked Gay-Berne units may also be used in the coarse-grained modelling 
of polymer chains including biological macromolecules 0, l^l • 

For illustrative purposes here, attention will focus on the "shape" term above, and the 
energy-dependent e terms will be taken to be constants. (Their contributions to forces and 
torques may be straightforwardly included using similar methods). It is therefore assumed 
that the potential may be written 



U = f/(r,a(f,a,b)) 



(8) 



The relation between the Gay-Berne distance function a and the geometry of ellipsoids has 



been clarified by Perram et al. 
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26 



(Xr, 



a 



. Define the matrices 



bm = '^bm®b„ 



(9) 



(but note that our A is A ^ in Ref. j27|). Here the and bm are semi-axis lengths of the 
ellipsoidal particles, and these have been used to define the un-normalised semi-axis vectors 



etc. These matrices may also be written in terms of the rotation matrices 



defined earher 



3 



mn 



Define the diagonal shape matrix S = diag(ai, 02, as), and the matrix a = Sa, whose rows 
are the un-normalised semi-axis vectors a^, i.e. [a]^^ = [dm]^ = o,mfi = cim.O'm.fj., m = 1,2, 3. 
This gives 

A = a"^S2a = a"^a . (11) 

Recall that the rotation matrix is orthogonal, a^a = aa^ = 1, whereas this is not true for a. 
The diameter a is defined through the following equations [3, ll^ 1^^ 

a = l/^, (12a) 
9? = £7-2 = if - (A + B)"^-f = if-H-i-f , (12b) 
where H = A + B . (12c) 

All of the a and b dependence lies within the matrix inverse while all the f dependence 
lies outside it. It is comparatively easy, though cumbersome, to evaluate the matrix inverse 
in a suitable form, perform the scalar products with r explicitly, and use eqns (El), @ for 

nnn ^ « 

the forces and torques [1^ Ha, |2^ . The following alternative approach was suggested by 
Perram et al. j^, and leads to more compact expressions based on eqns ((H), (j2I)- 



A. Scaled Potential 

First consider the special class of potentials for which 

U = U{rya'^{r, a, b)) = U{^) (13) 
where $ = r'^(f = r'^/a'^ = i r ■ (A + B) ~^ ■ r = | r ■ ■ r . (14) 

The simple form ((Tj) is an example. This form of potential is simple to implement, but 
it does have the disadvantage that the length scale is entirely determined by a; hence, a 
Lennard- Jones form will have an attractive range which is proportional to the repulsive 
diameter at any given orientation of molecules, which may be unphysical. The shifted form 
of the potential (see sec IIIIB|) was introduced by Gay and Berne [l^ to address this fault. 



On the other hand, in certain cases, e.g. for rapidly-varying, completely repulsive potentials, 



this may not 
Following 



2e_so much a problem. 

define an auxiliary vector k, by solving the equations 

H K = r ^ K = H^ r. (15) 



In what follows, there is no need to fully evaluate H ^. In terms of k 



\r ■ K , ^ = H^-r = K;. 

This last formula gives the interparticle force from 

dUd^ df/ 
^ = -^^ = -^''- ^'^^ 

To obtain the torque, write 

v-^ dU V- dU dU ^ (9$ 

= - ^ a„ X ^ = - )^ X ^ = - - ^ X ^ . (17) 

mm m 

The evaluation of d^/dam proceeds as follows. 

a$ , OH-' , OH , OA 
7) = 1'^--^ r = -ir-H H ■r = -|K-- k. 

The definition of A implies dA/dam^ = (g) + flm ® e^, where e^, fi = x,y, z is the 
appropriate unit vector in the space-fixed Cartesian coordinate system. Therefore 

d^/dam^, = -{k ■ am)K^, =^ / da.„, = -{n ■ arn)i^ ■ 

Combining with eqn ()17|) . and applying a similar formula to molecule B, gives the torques 

dU sr-^ dU 
ta = 2_^[i^ ■ am){am X K) = -^{i^ ■ A X k) , (18a) 

m 

tb = ^ ^(k; ■ b„,){brn X k) = ^(k • B X k) . (18b) 



These equations are essentially those of Perram et al. [22|. Using eqn ([THjl . it is possible to 
write them in the form 

= -K ■ A X / = rcA X / = {re - va) X f , (19a) 
tb = -A^ ■ B X / = rc-B X (-/) = {rc - r^) x (-/) (19b) 
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where rc is defined by 



rc = ta - 1^ - 1^ = rA + rcA (20a) 
= rB + B- K = rB + rcB (20b) 

(see Appendix ^ . The physical interpretation of this is that the intermolecular force may 
be taken to act at the "contact point" rc, and the torques about the molecular centres 
are then given by the usual moment formulae. Conservation of angular momentum follows 
immediately since ta + tb = (^ca - res) x / = -f ab >^ f = -r x f. 

It is illuminating to consider the case of identical uniaxial particles. Suppose the orien- 
tation of both particles is defined by unit vectors a^, along the 3-direction, and that the 
semiaxes are given by ai = 6i = 02 = &2 = Q-s = &3 = Then, the orthonormality 
relation Ylm ® ^"i = 1 allows the orientation matrices to be written in terms of as and 
63 alone: 

A = ill+ {q -il)a,®a,, B = f^l + (^J - £l)k ® k , 
H = A + B = 2ill + {q - il) (as ® as + bs ® ^3) . 

Inverting such a matrix is a standard exercise. One route is to observe that the orthonormal 
eigenvectors and eigenvalues are Q 

e± = , , A± = {il + il) ± {ij - il)a, . 6s , 



'2(l±as-bs) 
eo = X e_ , Aq = 2i1 . 

Then both H and its inverse may be expressed as sums: 

H = ^ Xm^m ® = AqI + (A+ - Ao)e+ (g) e+ + (A_ - Ao)e_ (g) e_ , 

m 

H"^ = J2 ^ra^m ® e„ = Aq 4 + (A;^ - Aq ^)e+ ® e+ + (Al^ - Aq ^)e_ ® e_ , 

m 

where the orthonormality property (g) = 1 has been used. 

It is convenient to define the elongation parameter x = (^y ~ ^i)/(^| + ^l); so that 
1 — X = 2£^/(£n + il); also scalar products of the relevant unit vectors: Ca = 0.3 ■ r, 
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Cfe = ^3 ■ r, Cab = as ■ 63; and to set Xc = XCab- Then 



H 



1 



2£l 
1 



X 



a3 (g) a3 + 63 ® 63 - Xc(a3 ® + 63 ® 0,3) 



X 



(a3 + bs) (g) {0,3 + 63) (03 - 63) O (03 - 63 



1 + Xc 



which gives the standard form '12'] for identical uniaxial particles 



a 



-2 _ 1 



f ■ H"^ ■ f = an'' <! 1 



(Ca + Cfe)2 _^ (Ca - CbY 



1 + Xc 



with a width parameter ctq = 2£_l. The auxiliary vector is 



K= 



■ r 



2il 



X 



1-x 



[Ca - XcCb)a3 + (Cfc - XcCa)b3 



The position vector of the "contact point" rc may be written 



rc = UrA + tb) - h 



X 



[Ca - XcCb)a3 - {Cb - XcCa)h 



' i-xl 

Having chosen the form of f/($), these last two equations make it easy to calculate forces 
and torques through eqns (|TI)|) and (|TT^ . 



B. Shifted Potential 

The shape- dependent part of the Gay-Berne potential Q| does not have the simple form 
of eqn (fT^. but uses the more general eqn with the definition of a given in eqn ()12b|) . 
This means that the potential is not a function of the single variable $ = jo^ ^ but instead 
depends on both r and a (or = cr~^) separately, typically through the shifted form 

Q = , (21) 

^min 

where cXmin is a constant. Define as before k = ■ r, = ■ ■ r = ^r^^r ■ k, so 

or 

This gives the interparticle force 

/ = - --^^^ ■ (1 - - ® = - r-'g^ r)r] • (22) 
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The torque derivation proceeds exactly as before, and the results are 

TA = r ■ A X k) , TB = r ■ B X k) . (23) 

Identifying part of the force = —r^^{dU/dip)K these equations may be written 

ta = -K ■ A X = rcA X , = ■ B X = res X (-/^) , (24) 

where Vc is given by eqn (j^UI) as before. Angular momentum conservation follows directly, 
since / — is parallel to f, so r x / = r x f^. 
As an example, consider the repulsive potential 

U = \ (25) 
[O, g'>2 

with g given by eqn pTj) and a given by eqn ()12b|) . It is sensible to set aram to be the 
minimum value of cr over all relative orientations. The potential then takes a constant value 
U = +eo when r a, for any orientation, and diverges as r a — cTmin- For example, 
for the identical uniaxial particles of the previous section, cimin = o"o = 2£_l for the prolate 
case i\i > i±, corresponding to side-by-side arrangements of the particles: but (Tmin = 2£|| for 
oblate, disk-hke, shapes, corresponding to the face-to-face orientation 29|. In the general 
case, eqn ()12b|) shows that if the minimum particle dimensions are amin = min(ai, 02, 03), 
bmin = min(6i, 62, 63), then 



'"mill 2i?inin 2(0^^;^ + b-^:^^) 



The geo,„etnca..y correct value for sphero.ds would be .„,., ^ a^-, but th,s d,sta,rce 

function does not faithfully represent the geometry of spheroids 27]. Then 

dU dUdgda , ,n o ,^ 

= -7r7r:r- = ^^^o 2^ - g ') ay2(x^i^ . 

dip og oa dip ^ ' 

Also 

or og or ^ 
From these two expressions, together with eqns (|22j) - (j24j) . the forces and torques may be 

easily calculated. 

The above equations generalise simply to the case of the full Gay-Berne potential, where 
the orientation-dependent energy functions ei and £2 are included [l7|. When the Lennard- 
Jones form ()25p is not truncated, the limiting long-range form is proportional to r~^, 

10 



as might be expected. However, a well-known criticism of the potential is that the full 
orient at ion- dependence is preserved at long range, whereas the van der Waals attractions 
in real molecules become isotropic at long range. It should be noted that Everaers and 
Ejtehadi js^ have derived a coarse-grained biaxial potential, called the RE-squared model, 
using Hamaker theory from colloid science. The result has features in common with the 
biaxial Gay-Berne potential, but also has the correct long-distance form. Babadi et al. 31 1 
have shown how to parametrise the RE-squared model using atomistic representations of 
small molecules. Forces and torques may be derived from this potential in a way analogous 
to that described above. 



IV. GAUSSIAN POTENTIAL 

Some years ago. Smith and Singer proposed an extension of the approach of Berne 
and Pechukas which they incorporated into the CCP5 Library program MDZOID. Each 
molecule is represented by a normalised Gaussian distribution of particle density, centred at 
va and r^, characterised by matrices A and B: 

= /o\iAi - ■ ■ - ^^)] ' (26a) 



where | ■ ■ ■ | is the determinant of a matrix. The interaction between volume elements dr^, 
dvb at Va and r^, respectively, is represented as a Gaussian function of their separation 



u 



[ra - n)dr„dn = Cr^ exp(-i |r, - nf /f) dr^dr, , (27) 



characterised by a range parameter i and a strength parameter C (with units of energy 
multiphed by volume). (More generally. Smith and Singer suggested fitting a desired 
atom- atom potential with a sum of such Gaussians, with positive or negative coefficients). 
The potential energy between molecules A and B is obtained by integrating u{ra — Tj,) over 
the density distributions PA{fa) and PB{fb)') the integrals may be performed analytically 
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(see Appendix IB)) and the result is, with r = r^B, 

U = C\H\-^/^ exp{-<^) , (28a) 
with H = £^1 + A + B , (28b) 
and $ = ir- (£2l + A + B)"^-r = ir--H-^-r. (28c) 

In the hmit £ — > 0, H and $ reduce to the definitions ()12c|) . (fT^ seen before. This potential 
belongs to the class of "core-softened" potentials, which have recently attracted attention 
in condensed matter physics. The interaction between polymer coils at low concentration, 
for instance, is well represented by a (spherically symmetric) Gaussian potential 33^ . 
The forces are derived exactly as for the scaled potential, sec. IIII Al 

with K = ■ r as before. The torques have contributions from the dependence of U on 
both $ and |H|. The $ contributions are given by 

= -K . A X / = qr^ X / , 
T% = -K-Bx f = q^x (-/) , 

which are eqns (fTIH) . with vca Qa — — A-k;, vcb — »• 9^ = B - k. The reason for the change 
in notation is that these equations do not define a contact point Vc as before in eqns (|^. 
because the new definition of H gives 

r + qA-qs = fi^, (30) 

rather than zero. Using the expression i|H|^^(9|H|/(9am) = H^^ ■ (see appendix IB]). the 
|H|-dependent terms in ta give 

^ipiE-^f^-IIHI-'f/E-^xllJI-f^E-xH-'-a™. (31) 

' ' m m m 

With the use of the Levi-Civita tensor e the final results may be expressed 

TA = U^amXH-^ -am + QA^ f = Ue: H ^ + qA^ f , (32a) 

m 

TB = Uj2bmy<H-' -bm-qBX f = U € : H'' B + Q ^ X {- f) . (32b) 
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Smith and Singer j32| obtained the same result, by a different method and written in a dif- 
ferent form (see Appendix^ll- Once more, angular momentum conservation follows straight- 
forwardly: 

TA + TB + r X f = Ue:H-\A + B) + {r + q^-qB) X f 

using eqn (jHOJ. The tensor 1 — £^H~^ is symmetric, so the double contraction with e gives 
zero, and f = Uk, so the last cross product also gives zero. 

V. ELLIPSOID CONTACT POTENTIAL 

As mentioned in section 1111 Bl the distance function appearing in the Gay-Berne and 
related potentials does not correctly represent the geometry of ellipsoidal particles. Formulae 



reflecting the true ellipsoid geometry were set out by Perram et al. 
takes the more general form 



The distance function 



= a-^ = a/3 r ■ {aA + (3B) ^ ■ r = af3 r ■ G'^ ■ r , (33) 
where a + (3 = 1. The matrix 

G = aA + /3B (34) 

plays a similar role to H in the previous sections; the case a = /5 = | corresponds to the 
standard Gay-Berne, or Gaussian overlap, definition of diameter seen earlier, and in this 
case G = ^H. The 'true' ellipsoid geometry corresponds to the case where a = (1 — /?) has 
been chosen to maximise ip, i.e. minimise a. This procedure is carried out numerically in 
the simulation. 

The "scaled potential" analogous to eqn (fT!?jl 

U = f/(rV(T2(f,a,b)) = f/($) , 
where $ = r'^ip = ja^ = a/Sr ■ (aA + (3B) ■ r = af3r ■ G^ ■ r , (35) 



is termed an ellipsoid contact potential j27|. The forces and torques are derived in a manner 
analogous to that of Section Fill Al Define an auxiliary vector A; by solving the equations 

G k = r ^ k = G^ r. (36) 
13 



In terms of k 



<^ = a(3r-k , (37) 
(9$ 

— = 2apG-^ ■ r = 2a(3k . (38) 
or 

This last formula will give the interparticle force from 

/ = -7^ = 2«/3fc . (39) 

The derivation of torques follows the pattern of previous sections. The results are 

TA = ^2a2/5^(fc . a^){am x fc) = ^2a^(3{k ■ A x fc) , (40a) 

m 

TB = ^2«/?2 ^(fc ■ K-,){brn X k) = ^2a(3\k ■ B X k) . (40b) 

m 

Once again these may be written 

TA = -ak ■ Ax f = rcA x / , (41a) 

= ■ B X / = rcB X (-/) , (41b) 

where Vc is now defined by 

rc = ta- ■ k = vb + ■ k (42) 

see Appendix E}. Conservation of angular momentum follows immediately. One final point 
27l |: the above derivatives take no account of the variation of a, /? as the particles are moved; 

but these extra terms vanish because the relevant functions have been minimised/maximised 

with respect to variations of these parameters. 

The corresponding formulae for more general (for example, shifted) potentials based on 

true ellipsoid geometry are 

(p = a(3r ■ ■ r = al3r~^r ■ k , (43) 

^ = 2af3G-^ ■ r = 2al3r-^k . (44) 
or 

The forces and torques are 

/ = -— r - — r-i-^ ■ l-r®r) = -—r - f[k-{k-r)r) , (45a 

Or Oip or or oip r^ 

TA = J2^k ■ am){am xk) = ^^^(fc ■ A X fc) , (45b) 

m 



0*09 ^ — ^ Oip 

m 



'^^ = ■ ^™)(^- ^) = a^^^^ ■ B X fc) . (45c) 

m 
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Identifying part of the force = —{dU/dip){2aP/r^)k the torques become 

TA = -ak . Ax fi^ = rcAX fk , (46a) 
TB = -f3k-Bx f, = rcB X (-/fc) , (46b) 

and, since f — f^. is parallel to f, angular momentum conservation follows as before. 

A modification of the ellipsoid contact potential, based on the surface- surface distance 
calculated along a vector parallel to the line of centres, has recently been proposed js^. 
Analytical formulae for the forces and torques are given in that paper. 



VI. CONCLUSIONS 



In the preceding sections it has been shown how compact expressions for forces and 
torques may be derived for pairwise intermolecular potentials which are written as functions 
of the molecular orientation matrices and the centre-centre vector. Illustrations have been 
given for various members of the "extended Gay-Berne family" ; in particular, for the Gaus- 
sian potential of Smith and Singer js^, much simpler expressions for the torques may be 
obtained in this way. 

The above expressions for forces and torques may be used in molecular dynamics sim- 
ulations, as well as hybrid Monte Carlo simulations or other applications. The y m ay be 
straightforwardly combined with symplectic and reversible integration algorithms based 
either on quaternion parameter representations of the rigid-body orientations, or on the 
rotation matrix itself 
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APPENDIX A: THE CONTACT FUNCTION 



This appendix recalls the relations leading to the definition of the "contact point" of 
eqn ()42|) . and the link with the Berne- Pechukas distance parameter. The following function 
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has a unique minimum as r varies for fixed a and (3 = 1 — a 25|, |26| : 

^(r; a, 13) = /3{r - va) ■ A"^ • (r - r^) + a{r - re) ■ ■ {r - Vb) ■ (Al) 



The location of this minimum, r = Vq, satisfies 
1 (9^ 



2 dr 



/3A-^ ■ {rc -rA) + aB ^ ■ {rc - r^) = . 



(A2) 



r=r0 



The definition G-fc = rAB = f, eqn (jHSl), with G = a;A+/3B, eqn (jHH), leads to an expression 
for rc 

r = G- fc = «A-fc + /3B-fc = (r^ - rc) + (^c - r^) , 
where = — aA ■ fc = + /3B ■ fc , 

which satisfies eqn f)A2|) . The value at the minimum is 

*lr-=rc = (^(^c - rA) ■ {-ak) + a{rc - Vb) ■ {(3k) 
= a(3r ■ k = a(3r ■ ■ r = $ , 

as defined in eqn ()37|) . Having identified rc and the function \& may be re- written 

^{r; a, P) = {r - rc) ■ {(3^-^ + aB'^) • (r - r^) + $ , (A3) 

as may be checked by expanding and comparing coefficients. 

The corresponding equations for the case a = /5 = | are conveniently written with the 



= rAB 


= r, eqn H = A 


+ B 


, eqn ^I^, 




^{r) = 


i(r-r-A) ■ A~^ ■ (r- 


rA) 


+ \{r-rB)-B-^ -{r-rB) , 


(A4) 


H K = 


A ■ K + B ■ K = {rA — 


rc) 


+ {rc - rB) , 


(A5) 


rc = 


- A ■ K = + B ■ 


K , 




(A6) 


\r=rc 


Ir ■ K = ^r ■ ■ r = 


E $ 




(A7) 




i(r-r-c)- (A-i + B- 




(r - rc) + $ . 


(A8) 



The generalised Berne-Pechukas overlap function is computed in terms of the two-centre 
Gaussian overlap integral 



I = J dr pA{r - rA)pB{r - rB) 



dr exp(— vl'(r)) 



exp(-$) 



Vs^v^aTb] 



(A9) 
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where the densities p^, Pb, are defined in eqns and the function \l/(r) is defined in 
eqn ()A4|) . The last step is carried out either by using Fourier transforms or by re- writing 
\l/(r) in the form of eqn ()A8|) and shifting origin to rc 27]; $ is defined in eqn ()A7|) . In 
practice, the Gaussian form of eqn ()A9|) is not commonly used as a potential in simulations; 
instead, the function $ is extracted from it, used to define a range parameter a, through 
(y9 = $/r^ and cr = eqns fll2j) of the text; and this range parameter is used in a variety 

of potentials. 

APPENDIX B: THE GAUSSIAN POTENTIAL 

The Gaussian potential of Smith and Singer is computed in terms of the double 



= y y dvadrb PA{ra - rA)pB{rb - rB)u{ra - Vb) , 

where u{ra — Vb) is defined by eqn (j2Zj). The two integrals are carried out successively: each 
has the same two-centre form as eqn ()A9|) . so they are straightforward. The result is 

U = C|H|-^/2 exp(-ir • H"^ ■ r) , where r = - and H = ^^1 + A + B , 

and this is eqn (j^Hj) of the main text. Smith and Singer derived an expression for the torque 
arising from the Gaussian potential, from first principles, treating uiva — Vb) as a site-site 
potential and integrating Vab x f^b over the density distributions. In the notation of sec IIV[ 
their result for the torque on A is 



where the definition qr^ = —A ■ n = — AH~ ■ r is used to simplify. Smith and Singer 
define the operation * to give a vector from two matrices (A * B)k = eij^AuBji (note the 
conventional summation over all three repeated indices i, j and /). This may be recognised 
as a double contraction of the Levi-Civita tensor e with a matrix product: 




integral 




ta = -{i^l + B) * (g^ ® r + ® qr^ + A - AH fC)U 
= - (H - A)"' * [qa ® (1 - AH-^) ■ r + (1 - AH-^)A] U , 



(Bl) 



A* B 



e : AB^ = e : BA^ . 



Using this, and noting that 



(H - A)-i(l - AH"^) = (H - A)-i(H - A)H 



H 



-1 
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eqn ()B1|) simplifies considerably: 



TA = e: (H"i ■ r ® + H-^A)f/ = e : ® + H-^A)f/ 



This is eqn (I32ap of the main text. 

In the derivation of eqn ()31|) . it is necessary to obtain the derivative of |H| with respect 
to a component of one of the rotation matrices. Use 



IHI' 



TJl ZJf Til 

XX xy xz 

Hyx Hyy Hy2^ 

Hzx ^zy Hzz 



+ 



Hxx Hxy Hxz 

TTI TJl TTI 

yx yy yz 

Hzx Hzy Hzz 



+ 



Hxx Hxy 

Hyx Hyy Hy^ 

H'zx H'^y H'^z 



Here, H^^ 



''5^u + Em « 



+ J^n^nfibnu, SO dHf,^/dc 



d\H\ 

dcirn.cf 





(^my 


^mz 




Hxx 


Hxy 


Hxz 




Hxx 


Hxy 


Hxz 


Hyx 


Hyy 


Hyz 


+ 


^my 








+ 


Hyx 


Hyy 


Hyz 


Hzx 


Hzy 


Hzz 




Hzx 


Hzy 


Hzz 




^mz 








'^^mx 


Hyy 


Hyz 




Hxy 


Hxz 




Hxy 


Hxz 




Hzy 


Hzz 






Hzy 


Hzz 






Hyy 


Hyz 



using the fact that H is symmetric. Moreover, each (signed) 2x2 determinant here is an 
element of the (signed) cofactor matrix, the transpose of which is multiplied by |H|. 
This allows the desired derivative to be written 



d\H\ 



2|H| ^OmsfH ^j^.^, + Omy[H ^] + «m2 [H ^] 



Collecting together results for all the components gives 



ar 



This is used in the derivation of eqn (jHH). 
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